1 IMPORT

1.1 Import librairies

library(ggplot2)
library(cowplot)
library(tidyverse)
library(dplyr)
library(Seurat)

setwd("~/Projects/HumanThymusProject/")
source("~/Projects/HumanThymusProject/scripts/colors_universal.R")

1.2 Import data

# import seurat objects
seur.nkt  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_inkt.rds")
seur.mait <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_mait.rds")
seur.gdt  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_gdt.rds")
seur.cd4  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd4.rds")
seur.cd8  <- readRDS("./data_github/seurat_objects_thymus/seurat_thymus_cd8.rds")

# sanity check seurat objects correspond to Fig 2
DimPlot(seur.nkt,  group.by="clusters_per_lineage", label=T, cols=cols_thym_nkt)

DimPlot(seur.mait, group.by="clusters_per_lineage", label=T, cols=cols_thym_mait)

DimPlot(seur.gdt,  group.by="clusters_per_lineage", label=T, cols=cols_thym_gdt)

# import GEP usage dataframe
gep_usage <- read.table("./data_github/cNMF/non_imputed_cNMF_allcells.usages.k_12.dt_0_02.consensus.txt", header=T)
colnames(gep_usage) <- paste0("GEP", c(2,5,3,1,4,12,6,7,8,10,9,11), "_usage")
gep_usage <- gep_usage[,paste0("GEP", 1:12, "_usage")]
head(gep_usage)
##                     GEP1_usage GEP2_usage GEP3_usage GEP4_usage  GEP5_usage
## CD4_1_Thymus_878524  0.0000000 0.00000000 0.01697069 0.02701780 0.003293012
## CD4_1_Thymus_679405  0.0000000 0.00000000 0.00000000 0.00000000 0.000000000
## CD4_1_Thymus_281188  0.0000000 0.01499587 0.02538062 0.01897941 0.000000000
## CD4_1_Thymus_285601  0.3893558 0.19701648 0.01374016 0.03301096 0.009603295
## CD4_1_Thymus_174537  0.0000000 0.00000000 0.01813514 0.02750585 0.001959355
## CD4_1_Thymus_128406  0.0000000 0.02029776 0.02094953 0.01324090 0.001427132
##                     GEP6_usage   GEP7_usage GEP8_usage GEP9_usage GEP10_usage
## CD4_1_Thymus_878524          0 0.0000000000 0.01670451 0.06574529  0.14537527
## CD4_1_Thymus_679405          0 0.0000000000 0.00000000 0.00000000  0.00000000
## CD4_1_Thymus_281188          0 0.0002823688 0.01075830 0.09749147  0.21138272
## CD4_1_Thymus_285601          0 0.0178349300 0.22377394 0.05961335  0.01831705
## CD4_1_Thymus_174537          0 0.0000000000 0.01610514 0.06646363  0.15172973
## CD4_1_Thymus_128406          0 0.0055636405 0.01326269 0.08876818  0.21684438
##                     GEP11_usage GEP12_usage
## CD4_1_Thymus_878524 0.655783300  0.00000000
## CD4_1_Thymus_679405 1.147190900  0.00000000
## CD4_1_Thymus_281188 0.563456800  0.00000000
## CD4_1_Thymus_285601 0.008145613  0.03549316
## CD4_1_Thymus_174537 0.649509700  0.00000000
## CD4_1_Thymus_128406 0.565640400  0.00000000

2 ANALYSIS

2.1 Add GEP1-GEP11 usage to seurat objects

# check that rownames are in same order
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data),])==rownames(seur.cd4@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data),])==rownames(seur.cd8@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data),])==rownames(seur.nkt@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data),])==rownames(seur.mait@meta.data), useNA="ifany")
# table(rownames(gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data),])==rownames(seur.gdt@meta.data), useNA="ifany") # all TRUE

# add GEP usage to seurat objects
gep_colnames <- paste0("GEP", 1:11, "_usage")
seur.cd4@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd4@meta.data), gep_colnames]
seur.cd8@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.cd8@meta.data), gep_colnames]
seur.nkt@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.nkt@meta.data), gep_colnames]
seur.mait@meta.data[,gep_colnames] <- gep_usage[rownames(gep_usage) %in% rownames(seur.mait@meta.data), gep_colnames]
seur.gdt@meta.data[,gep_colnames]  <- gep_usage[rownames(gep_usage) %in% rownames(seur.gdt@meta.data), gep_colnames]

2.2 Plot GEP usage per thymic lineage

# Plot
seur_vector <- list("iNKT"=seur.nkt, "MAIT"=seur.mait, "GDT"=seur.gdt, "CD4"=seur.cd4, "CD8"=seur.cd8)

plist_massive <- list()
for(gep in gep_colnames){
  # print(gep)
  # get row of plots (for one GEP)
  plist <- list()
  for(i_seur in names(seur_vector)){
    # print(i_seur)
    seur <- seur_vector[[i_seur]]
    # get legend once
    if(gep=="GEP1_usage" & i_seur=="iNKT"){
      plegend <- ggpubr::get_legend(
        SCpubr::do_FeaturePlot(seur,  features=gep, ncol=1, legend.position="right", legend.title="GEP usage", border.size=1, pt.size=2, order=T)+
          scale_color_viridis_c(limits=c(0,1.25), option="B")+
          theme(plot.background = element_rect(color = "black"))
      )
    }
    # make plot
    legendpos <- "none"
    plt_title <- ""
    if(gep=="GEP1_usage"){plt_title=i_seur}
    p <- SCpubr::do_FeaturePlot(seur,  features=gep, ncol=1, legend.position=legendpos, border.size=1, pt.size=2, order=T, plot.title=plt_title)+
      scale_color_viridis_c(limits=c(0,1.25), option="B")+
      theme(plot.background = element_rect(color = "black", linewidth = 2),
            plot.margin=unit(c(0,5.5,0,5.5), "points"))
    plist[[i_seur]] <- ggrastr::rasterise(p, layers="Point", dpi=300)
  }
  prow_gep <- plot_grid(plotlist=plist, nrow=1, scale=0.9)
  plist_massive[[gep]] <- prow_gep 
}

plot_grid(
  plot_grid(plotlist=plist_massive, ncol=1, align="h"),
  ggpubr::as_ggplot(plegend), ncol=2, rel_widths = c(5, .5), scale=c(0.95, 5)
)

3 SESSION INFO

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] SeuratObject_4.1.3 Seurat_4.3.0.1     forcats_1.0.0      stringr_1.5.0     
##  [5] dplyr_1.1.2        purrr_1.0.1        readr_2.1.4        tidyr_1.3.0       
##  [9] tibble_3.2.1       tidyverse_1.3.2    cowplot_1.1.1      ggplot2_3.4.2     
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.4.2           backports_1.4.1        plyr_1.8.8            
##   [4] igraph_1.5.0           lazyeval_0.2.2         sp_2.0-0              
##   [7] splines_4.1.3          listenv_0.9.0          scattermore_1.2       
##  [10] digest_0.6.32          htmltools_0.5.5        viridis_0.6.3         
##  [13] fansi_1.0.4            magrittr_2.0.3         tensor_1.5            
##  [16] googlesheets4_1.1.1    cluster_2.1.4          ROCR_1.0-11           
##  [19] tzdb_0.4.0             globals_0.16.2         modelr_0.1.11         
##  [22] matrixStats_1.0.0      timechange_0.2.0       spatstat.sparse_3.0-2 
##  [25] colorspace_2.1-0       rvest_1.0.3            ggrepel_0.9.3         
##  [28] haven_2.5.3            xfun_0.39              crayon_1.5.2          
##  [31] jsonlite_1.8.7         progressr_0.13.0       spatstat.data_3.0-1   
##  [34] survival_3.5-5         zoo_1.8-12             glue_1.6.2            
##  [37] polyclip_1.10-4        gtable_0.3.3           gargle_1.5.1          
##  [40] leiden_0.4.3           car_3.1-2              future.apply_1.11.0   
##  [43] abind_1.4-5            scales_1.2.1           DBI_1.1.3             
##  [46] SCpubr_1.1.2           rstatix_0.7.2          spatstat.random_3.1-5 
##  [49] miniUI_0.1.1.1         Rcpp_1.0.10            viridisLite_0.4.2     
##  [52] xtable_1.8-4           reticulate_1.30        htmlwidgets_1.6.2     
##  [55] httr_1.4.6             RColorBrewer_1.1-3     ellipsis_0.3.2        
##  [58] ica_1.0-3              pkgconfig_2.0.3        farver_2.1.1          
##  [61] sass_0.4.6             uwot_0.1.16            dbplyr_2.3.2          
##  [64] deldir_1.0-9           utf8_1.2.3             tidyselect_1.2.0      
##  [67] labeling_0.4.2         rlang_1.1.1            reshape2_1.4.4        
##  [70] later_1.3.1            munsell_0.5.0          cellranger_1.1.0      
##  [73] tools_4.1.3            cachem_1.0.8           cli_3.6.1             
##  [76] generics_0.1.3         broom_1.0.5            ggridges_0.5.4        
##  [79] evaluate_0.21          fastmap_1.1.1          yaml_2.3.7            
##  [82] goftest_1.2-3          knitr_1.43             fs_1.6.2              
##  [85] fitdistrplus_1.1-11    RANN_2.6.1             pbapply_1.7-2         
##  [88] future_1.33.0          nlme_3.1-162           mime_0.12             
##  [91] ggrastr_1.0.2          xml2_1.3.4             compiler_4.1.3        
##  [94] rstudioapi_0.14        beeswarm_0.4.0         plotly_4.10.2         
##  [97] png_0.1-8              ggsignif_0.6.4         spatstat.utils_3.0-3  
## [100] reprex_2.0.2           bslib_0.5.0            stringi_1.7.12        
## [103] highr_0.10             lattice_0.21-8         Matrix_1.5-4.1        
## [106] vctrs_0.6.3            pillar_1.9.0           lifecycle_1.0.3       
## [109] spatstat.geom_3.2-1    lmtest_0.9-40          jquerylib_0.1.4       
## [112] RcppAnnoy_0.0.21       data.table_1.14.8      irlba_2.3.5.1         
## [115] httpuv_1.6.11          patchwork_1.1.2        R6_2.5.1              
## [118] promises_1.2.0.1       KernSmooth_2.23-21     gridExtra_2.3         
## [121] vipor_0.4.5            parallelly_1.36.0      codetools_0.2-19      
## [124] assertthat_0.2.1       MASS_7.3-60            withr_2.5.0           
## [127] sctransform_0.3.5      parallel_4.1.3         hms_1.1.3             
## [130] grid_4.1.3             rmarkdown_2.23         carData_3.0-5         
## [133] googledrive_2.1.1      Cairo_1.6-0            Rtsne_0.16            
## [136] ggpubr_0.6.0           spatstat.explore_3.2-1 shiny_1.7.4           
## [139] lubridate_1.9.2        ggbeeswarm_0.7.2